(in-package :differentiator)

(defun diffn (expr dx)
  "Differentiate EXPR with respect to DX. Make no attempt to
   perform constant folding or other simplification. EXPR must be
   in normalized form: all operations are either unary or
   binary."
  (etypecase expr
    (number 0)
    (symbol (if (eql expr dx) 1 0))
    (list
     (destructuring-bind (op . args) expr
       (ecase op
	 (sin `(* (cos ,@args) ,(diffn (car args) dx)))
	 (cos `(* (* -1 (sin ,@args)) ,(diffn (car args) dx)))
	 (+ `(+ ,(diffn (first args) dx) ,(diffn (second args) dx)))
	 (- `(- ,(diffn (first args) dx) ,(diffn (second args) dx)))
	 (* `(+ (* ,(first args) ,(diffn (second args) dx))
		(* ,(diffn (first args) dx) ,(second args))))
	 (/ `(- (/ ,(diffn (first args) dx) ,(second args))
		(/ (* ,(first args) ,(diffn (second args) dx))
		   (expt ,(second args) 2))))
	 (expt `(* (* ,(second args) (expt ,(first args) ,(1- (second args))))
		   ,(diffn (first args) dx)))
	 (exp `(exp ,@args)))))))

(defun normalize (expr)
  (etypecase expr
    ((or number symbol) expr)    
    (list
     (destructuring-bind (op . args) expr
       (let ((nargs (mapcar #'normalize args)))
	 (if (< 2 (length nargs))
	     (if (or (eql op '+)
		     (eql op '*))
		 (reduce #'(lambda (a b) (list op a b)) nargs :from-end t)
		 (list op
		       (car nargs)
		       (reduce #'(lambda (a b) (list (if (eql op '-) '+ '*) a b)) (cdr nargs) :from-end t)))
	     (cons op nargs)))))))

(defun diff (expr dx)
  "Differentiate EXPR with respect to DX. Make no attempt to
   perform constant folding or other simplification."
  (diffn (normalize expr) dx))